﻿#region netDxf library licensed under the MIT License
// 
//                       netDxf library
// Copyright (c) Daniel Carvajal (haplokuon@gmail.com)
// 
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
// 
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
// 
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
// 
#endregion

// This is a translation to C# from the original C++ code of the Geometric Tool Library
// Original license
// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2022
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 6.0.2022.01.06

using System;
using System.Linq;

namespace netDxf.GTE
{
    public static class Integration
    {
        // A simple algorithm, but slow to converge as the number of samples
        // is increased.  The 'numSamples' needs to be two or larger.
        public static double TrapezoidRule(int numSamples, double a, double b, Func<double, double> integrand)
        {
            double h = (b - a) / (numSamples - 1.0);
            double result = 0.5 * (integrand(a) + integrand(b));
            for (int i = 1; i <= numSamples - 2; ++i)
            {
                result += integrand(a + i * h);
            }
            result *= h;

            return result;
        }

        // The trapezoid rule is used to generate initial estimates, but then
        // Richardson extrapolation is used to improve the estimates.  This is
        // preferred over TrapezoidRule.  The 'order' must be positive.
        public static double Romberg(int order, double a, double b, Func<double, double> integrand)
        {
            double half = 0.5;
            double[][] rom = new double[order][];
            for (int i = 0; i < order; i++)
            {
                rom[i] = new double[2];
            }

            double h = b - a;
            rom[0][0] = half * h * (integrand(a) + integrand(b));
            for (int i0 = 2, p0 = 1; i0 <= order; i0++, p0 *= 2, h *= half)
            {
                // Approximations via the trapezoid rule.
                double sum = 0.0;
                int i1;
                for (i1 = 1; i1 <= p0; ++i1)
                {
                    sum += integrand(a + h * (i1 - half));
                }

                // Richardson extrapolation.
                rom[0][1] = half * (rom[0][0] + h * sum);
                for (int i2 = 1, i2m1 = 0, p2 = 4; i2 < i0; i2++, i2m1++, p2 *= 4)
                {
                    rom[i2][1] = (p2 * rom[i2m1][1] - rom[i2m1][0]) / (p2 - 1);
                }

                for (i1 = 0; i1 < i0; ++i1)
                {
                    rom[i1][0] = rom[i1][1];
                }
            }

            return rom[order - 1][0];
        }

        private readonly struct Info
        {
            public readonly int NumBits;
            public readonly double[] Product;

            public Info(int numBits, double[] product)
            {
                this.NumBits = numBits;
                this.Product = product;
            }
        };

        // Gaussian quadrature estimates the integral of a function f(x)
        // defined on [-1,1] using
        //   integral_{-1}^{1} f(t) dt = sum_{i=0}^{n-1} c[i]*f(r[i])
        // where r[i] are the roots to the Legendre polynomial p(t) of degree
        // n and
        //   c[i] = integral_{-1}^{1} prod_{j=0,j!=i} (t-r[j]/(r[i]-r[j]) dt
        // To integrate over [a,b], a transformation to [-1,1] is applied
        // internally:  x - ((b-a)*t + (b+a))/2.  The Legendre polynomials are
        // generated by
        //   P[0](x) = 1, P[1](x) = x,
        //   P[k](x) = ((2*k-1)*x*P[k-1](x) - (k-1)*P[k-2](x))/k, k >= 2
        // Implementing the polynomial generation is simple, and computing the
        // roots requires a numerical method for finding polynomial roots.
        // The challenging task is to develop an efficient algorithm for
        // computing the coefficients c[i] for a specified degree.  The
        // 'degree' must be two or larger.
        public static void ComputeQuadratureInfo(int degree, out double[] roots, out double[] coefficients)
        {
            const double zero = 0.0;
            const double one = 1.0;
            const double half = 0.5;
            double[][] poly = new double[degree + 1][];

            poly[0] = new[] {zero};

            poly[1] = new[] {zero, one};

            for (int nm = 2, nm1 = 1, nm2 = 0, np1 = 3; nm <= degree; nm++, nm1++, nm2++, np1++)
            {
                double mult0 = nm1 / (double) nm;
                double mult1 = (2.0 * nm - 1.0) / nm;

                poly[nm] = new double[np1];
                poly[nm][0] = -mult0 * poly[nm2][0];
                for (int i = 1, im1 = 0; i <= nm2; i++, im1++)
                {
                    poly[nm][i] = mult1 * poly[nm1][im1] - mult0 * poly[nm2][i];
                }
                poly[nm][nm1] = mult1 * poly[nm1][nm2];
                poly[nm][nm] = mult1 * poly[nm1][nm1];
            }

            RootsPolynomial.Find(degree, poly[degree], 2048, out roots);

            coefficients = new double[roots.Length];
            int n = roots.Length - 1;
            double[] subroots = new double[n];
            for (int i = 0; i < roots.Length; i++)
            {
                double denominator = 1.0;
                for (int j = 0, k = 0; j < roots.Length; j++)
                {
                    if (j != i)
                    {
                        subroots[k++] = roots[j];
                        denominator *= roots[i] - roots[j];
                    }
                }

                double[] delta = {-one - subroots.Last(), one - subroots.Last()};
                double[][] weights = new double[n][];
                weights[0] = new[] {half * delta[0] * delta[0], half * delta[1] * delta[1]};

                for (int k = 1; k < n; k++)
                {
                    double dk = k;
                    double mult = -dk / (dk + 2.0);
                    weights[k][0] = mult * delta[0] * weights[k - 1][0];
                    weights[k][1] = mult * delta[1] * weights[k - 1][1];
                }

                int numElements = 1 << (n - 1);
                Info[] info = new Info[numElements];
                info[0] = new Info(
                    0,
                    new[]
                    {
                        one,
                        one
                    }
                );

                for (int ipow = 1, r = 0; ipow < numElements; ipow <<= 1, r++)
                {
                    info[ipow] = new Info(
                        1,
                        new[]
                        {
                            -one - subroots[r],
                            +one - subroots[r]
                        }
                    );

                    for (int m = 1, j = ipow + 1; m < ipow; m++, j++)
                    {
                        info[j] = new Info(
                            info[m].NumBits + 1,
                            new[]
                            {
                                info[ipow].Product[0] * info[m].Product[0],
                                info[ipow].Product[1] * info[m].Product[1]
                            }
                        );
                    }
                }

                double[][] sum = new double[n][];
                for (int ni = 0; ni < n; ni++)
                {
                    sum[ni] = new[] {zero, zero};
                }

                for (int k = 0; k < info.Length; ++k)
                {
                    sum[info[k].NumBits][0] += info[k].Product[0];
                    sum[info[k].NumBits][1] += info[k].Product[1];
                }

                double[] total = {zero, zero};
                for (int k = 0; k < n; ++k)
                {
                    total[0] += weights[n - 1 - k][0] * sum[k][0];
                    total[1] += weights[n - 1 - k][1] * sum[k][1];
                }

                coefficients[i] = (total[1] - total[0]) / denominator;
            }
        }

        public static double GaussianQuadrature(double[] roots, double[] coefficients, double a, double b, Func<double, double> integrand)
        {
            const double half = 0.5;
            double radius = half * (b - a);
            double center = half * (b + a);
            double result = 0.0;
            for (int i = 0; i < roots.Length; i++)
            {
                result += coefficients[i] * integrand(radius * roots[i] + center);
            }
            result *= radius;
            return result;
        }
    };
}

